rm(list=ls())

library(foreign)
library(car)
library(stargazer)
library(lattice)
library(ggplot2)
library(texreg)
library(psych)
library(sjPlot)
library(sjmisc)
library(ggplot2)
library("gridExtra")
library(lme4)



setwd("/Users/amandaclayton/Dropbox/Party discipline/Draft/Submissions/APSR/R&R!/R&R 2/Accepted!!/Replication Files")


#read data
modeldata<-read.csv("modeldata.csv", header=T, stringsAsFactors=FALSE, fileEncoding="latin1")


#number of women in each country 

sort(tapply(modeldata$Female, modeldata$Country, sum))

##Table 2: Descriptive Statistics 

t.test(scores ~ Female, data=modeldata)
t.test(Ruling ~ Female, data=modeldata)
t.test(MPYears ~ Female, data=modeldata)
t.test(Minister ~ Female, data=modeldata)
t.test(PartyLeader ~ Female, data=modeldata)
t.test(Edu ~ Female, data=modeldata)
t.test(PoliticalJob ~ Female, data=modeldata)


##Table 3, main regression results, DV is discipline scores  

#Note to self, update table and text 

model1<-lm(scores~Female + Party, data=modeldata)
summary(model1)  
screenreg(model1)

model2<-lm(scores~Female + Ruling + MPYears+ Minister+ PartyLeader + Edu + PoliticalJob + Party, data=modeldata)
summary(model2) 

screenreg(list(model1, model2), stars = c(0.01, 0.05, 0.1), digits = 3)

texreg(list(model1, model2), stars = c(0.01, 0.05, 0.1), digits = 3)


#Figure 1 

#load speech data
hansard<-read.csv("HansardData.csv")

#basic gender differences referenced in text

lm1<-lm(numb_stand ~ Female + party, data = hansard)
summary(lm1)

lm2<-lm(self_ref_perc ~ Female + party, data = hansard)  
summary(lm2)

#As noted in text: analysis excludes Zimbabwe b/c quota imposed subsequent year, higher than usual re-election rates for women incumbents

sub<-subset(hansard, Country!="Zimbabwe")

lm3<-lm(Reelected ~ Female + numb_stand + Female * numb_stand + party, data= sub)
summary(lm3)


##plot interaction of women over time for men and women 

set_theme(base = theme_grey())

set_theme(
  legend.inside = TRUE,         # legend inside plot
  legend.pos = "top right",  # legend position inside plot
 # legend.title.size = .8,
 # geom.label.size = 3
)


# make categorical
sub $temp<-as.factor(sub $Female)
sub $temp <-recode(sub $temp, "0 = 'Men' ; 1 = 'Women' ")
sub $temp <-to_factor(sub $temp)


fit1<-glm(Reelected ~temp+ numb_stand + temp * numb_stand +  party, data= sub)
summary(fit1)

plot_model(fit1, type = "pred", terms = c("numb_stand", "temp"), title = "", axis.title = c("Number of words (standardized)", "Reelected probability" ), legend.title = "MP Gender", colors="bw")

#individual correlations 

men<-subset(sub, Female==0)
women<-subset(sub, Female==1)

cor.test(men$numb_stand, men$Reelected, data=merge_sub) #positively correlated, sig at 0.003
cor.test(women $numb_stand, women $Reelected, data=merge_sub) #negative, not significat  


#Figure 2

#Note: Models correspond with Table A16 in SI 12 

model1<-lm(WomenRightsAgg~Female+scores+ MPYears+ Minister+ PartyLeader + Ruling+ Edu + PoliticalJob + Party + Country, data=modeldata)
summary(model1)

model2<-lm(WomenRightsAgg~Female+scores+MPYears+ Minister+ PartyLeader + Ruling+I(scores*Female)+ Edu + PoliticalJob + Party + Country, data=modeldata)
summary(model2)

screenreg(list(model1, model2))
texreg(list(model1, model2), stars = c(0.01, 0.05, 0.1), digits = 3)


set_theme(base = theme_grey())

set_theme(
  legend.inside = TRUE,         # legend inside plot
  legend.pos = "top right",  # legend position inside plot
 # legend.title.size = .8,
 # geom.label.size = 3
)

# make categorical
modeldata$Gender<-as.factor(modeldata$Female)
modeldata$Gender <-recode(modeldata$Gender, "0 = 'Men' ; 1 = 'Women' ")
modeldata$Gender <-to_factor(modeldata$Gender)

fit<-lm(WomenRightsAgg~Gender+scores+MPYears+ Minister+ PartyLeader + Ruling+scores*Gender+ Edu + PoliticalJob + Party + Country, data=modeldata)
summary(fit)

plot_model(fit, type = "pred", terms = c("scores", "Gender"), title = "", axis.title = c("Discipline Scores", "Women's Rights a Top Priority" ), legend.title = "MP Gender", colors="bw")


########################
##Key Appendix Tables and Images ##
########################

#discipline scores correlation with other indicators 

#Table A4
modeldata$logDM<-log(modeldata$DistMag)
correlation.matrix <- cor(modeldata[,c("scores","Proportional", "logDM", "Polity")]) 
stargazer(correlation.matrix, title="Correlation Matrix")

cor.test(modeldata$Proportional, modeldata$scores) 
cor.test(modeldata$logDM, modeldata$scores) 
cor.test(modeldata$Polity, modeldata$scores) 


##Table A5 in SI 6, interaction terms with previous political experience

model6<-lm(scores~Female + Ruling + I(Female * PoliticalJob) + MPYears+ Minister+ PartyLeader + Edu + PoliticalJob + Party, data=modeldata)
summary(model6)  

texreg(list(model6), stars = c(0.01, 0.05, 0.1), digits = 3)


##Table A7 and Table A8 SI, 7

lm1<-lm(numb_stand ~ Female + party, data = hansard)
summary(lm1)

texreg(list(lm1), stars = c(0.01, 0.05, 0.1), digits = 3)

lm2<-lm(self_ref_perc ~ Female + party, data = hansard)  
summary(lm2)

texreg(list(lm1, lm2), stars = c(0.01, 0.05, 0.1), digits = 6)


#Table A15, SI 11: Robustness, remove cases where women are less than 10 percent of country sample 

#number of women in each counry sample
sort(tapply(modeldata$Female, modeldata$Country, sum)) #

subsamp<-subset(modeldata, modeldata$Country != "Botswana" & modeldata$Country != "Ghana" & modeldata$Country != "Kenya" & modeldata$Country != "Benin" & modeldata$Country != "Burkina" & modeldata$Country != "Mali" )

model1<-lm(scores~Female + Party, data=subsamp)
summary(model1)  
screenreg(model1)

model2<-lm(scores~Female + Ruling + MPYears+ Minister+ PartyLeader + Edu + PoliticalJob + Party, data= subsamp)
summary(model2) #women sig less 
screenreg(list(model1, model2))

texreg(list(model1, model2), stars = c(0.01, 0.05, 0.1), digits = 3)


#Table A13 and A14, SI 10: Robustness, matrilineal v. patrilineal women 

t.test(scores~Female, data=subset(modeldata, modeldata$Matrilin ==1)) #still sig diff
t.test(scores~Female, data=subset(modeldata, modeldata$Matrilin ==0)) #still sig diff

lm1<-lm(scores ~ Female + Matrilin + I(Female* Matrilin), data=modeldata)
summary(lm1)


#SI 10: Cross-level interactions, Table A10 

mixed1<-lmer(scores ~ Female + Ruling+MPYears+ Minister+ PartyLeader + Edu + PoliticalJob + Quota +  Polity + log(GDPPC) + log(DistMag) + ABWomen + WomenRep + PREG  + (1 | Country), data=modeldata) #
summary(mixed1)

#interactions: 
mixed2<-lmer(scores ~ Female + Ruling+MPYears+ Minister+ PartyLeader + Edu + PoliticalJob + Quota +  Polity + log(GDPPC) + log(DistMag) + ABWomen + WomenRep + PREG + (Quota*Female) + (1 + Quota| Country), data=modeldata) # random intercept, cross-level interaction
summary(mixed2)

mixed3<-lmer(scores ~ Female + Ruling+MPYears+ Minister+ PartyLeader + Edu + PoliticalJob + Quota +  Polity + log(GDPPC) + log(DistMag) + ABWomen  + WomenRep + PREG + (log(DistMag)*Female) + (1 + log(DistMag)| Country), data=modeldata) # random intercept, cross-level interaction
summary(mixed3)

mixed4<-lmer(scores ~ Female + Ruling+MPYears+ Minister+ PartyLeader + Edu + PoliticalJob + Quota +  Polity + log(GDPPC) + log(DistMag) + ABWomen  + WomenRep +PREG + (Polity *Female) + (1 + Polity| Country), data=modeldata) # random intercept, cross-level interaction
summary(mixed4)

mixed5<-lmer(scores ~ Female + Ruling+MPYears+ Minister+ PartyLeader + Edu + PoliticalJob + Quota +  Polity + log(GDPPC) + log(DistMag) + ABWomen  + WomenRep +PREG + (WomenRep *Female) + (1 + WomenRep | Country), data=modeldata) # random intercept, cross-level interaction
summary(mixed5)

mixed6<-lmer(scores ~ Female + Ruling+MPYears+ Minister+ PartyLeader + Edu + PoliticalJob + Quota +  Polity + log(GDPPC) + log(DistMag) + ABWomen  + WomenRep +PREG + (ABWomen *Female) + (1 + ABWomen | Country), data=modeldata) # random intercept, cross-level interaction
summary(mixed6)

mixed7<-lmer(scores ~ Female + Ruling+MPYears+ Minister+ PartyLeader + Edu + PoliticalJob + Quota +  Polity + log(GDPPC) + log(DistMag) + ABWomen  + WomenRep + PREG + (PREG *Female) + (1 + PREG | Country), data=modeldata) # random intercept, cross-level interaction
summary(mixed7)

screenreg(list(mixed1, mixed2, mixed3, mixed4, mixed5, mixed6, mixed7), stars = c(0.01, 0.05, 0.1), digits = 3)

texreg(list(mixed1, mixed2, mixed3, mixed4, mixed5, mixed6, mixed7), stars = c(0.01, 0.05, 0.1), digits = 3)

#Changes in women's rights over time and in party discipline over time, among women

##SI Table A17, Figure A4

model3<-lm(WomenRightsAgg~ MPYears+ Minister+ PartyLeader + Ruling+ Edu + PoliticalJob + Party, data=modeldata)
summary(model3)

model4<-lm(scores~ MPYears+ Minister+ PartyLeader + Ruling+ Edu + PoliticalJob + Party, data= modeldata)
summary(model4)

screenreg(list(model3, model4))
texreg(list(model3, model4), stars = c(0.01, 0.05, 0.1), digits = 3)

##plot interaction of women over time for men and women 

fit1<-lm(WomenRightsAgg~Gender+MPYears+ Minister+ PartyLeader + Ruling+MPYears*Gender+ Edu + PoliticalJob + scores + Party, data=modeldata)
summary(fit1)

plot1 <-plot_model(fit1, type = "pred", terms = c("MPYears", "Gender"), title = "", axis.title = c("MP Years in Office", "Women's Rights a Top Priority" ), legend.title = "MP Gender", colors="bw")

#for party discipline 

fit2<-lm(scores~Gender+MPYears+ Minister+ PartyLeader + Ruling+MPYears*Gender+ Edu + PoliticalJob + Party, data=modeldata)
summary(fit2)

plot2 <- plot_model(fit2, type = "pred", terms = c("MPYears", "Gender"), title = "", axis.title = c("MP Years in Office", "Party Discipline Score" ), legend.title = "MP Gender", colors="bw")

grid.arrange(plot1, plot2, ncol=2)

